﻿#pragma once
#include <vector>
#include <map>
#include <memory>
#include <functional>
#include <iostream>
#include <algorithm>
using namespace std;
class Vector : public vector<double>
{
public:
	Vector(int n = 0){ (*this).resize(n); }
	Vector(int n, istream& f) //Ввод передаем имя файла, формат n v(0)...v(n) - размерность, далее элементы
	{
		(*this).resize(n);
		for (int i = 0; i < n; i++)
		{
			f >> (*this)[i];
		}
		//if (!f)throw exception("End of file read");
	}
	double operator*(const Vector &vec)const
	{
		double sum = 0;
		for (unsigned int i = 0; i < size(); i++)sum += (*this)[i] * vec[i];
		return sum;
	}
	void operator*=(double n){ for (auto& i : (*this))i *= n; }
	void operator+=(const Vector &vec){ for (unsigned int i = 0; i < size(); i++)(*this)[i] += vec[i]; }
	void operator-=(const Vector &vec){ for (unsigned int i = 0; i < size(); i++)(*this)[i] -= vec[i]; }
	Vector operator-(Vector &vec)const
	{
		Vector result;
		result.resize(vec.size());
		for (unsigned int i = 0; i < size(); i++)result[i] = (*this)[i] - vec[i];
		return result;
	}
	Vector operator+(Vector &vec)const
	{
		Vector result;
		result.resize(vec.size());
		for (unsigned int i = 0; i < size(); i++)result[i] = (*this)[i] + vec[i];
		return result;
	}
	Vector operator*(double n)const
	{
		Vector result;
		result.resize((*this).size());
		for (unsigned int i = 0; i < size(); i++)result[i] = (*this)[i] * n;
		return result;
	}

	inline friend static ostream& operator << (ostream& out, Vector& v)
	{
		for (auto i : v)out << i << " ";
		out << endl;
		return out;
	}
};

class MappedRowSkyLineFormat;
typedef MappedRowSkyLineFormat SpecialInputClass;

class Matrix
{
public:
	static shared_ptr<Matrix> CreateInstance(shared_ptr<SpecialInputClass> matr){};
	virtual void Run(function<void(int, int, double)>) const = 0;
	virtual void MultMatrixVector(const Vector& vec, Vector& res) const
	{
		res.resize(Size());
		for (auto& r : res)r = 0;
		Run([&](int i, int j, double d){res[i] += vec[j] * d; });
	}
	virtual void MultTransMatrixVector(const Vector& vec, Vector& res) const = 0;
	virtual void multLTransInvToVector(bool withDiag, const Vector& in, Vector& out) const// (L [+D])^T)^-1^ * vector, нужно для бисопряженных градиентов
	{
		throw exception("У предобуславливателя хочу чтобы было multLTransInvToVector, а она есть только в одном классе");
	}

	virtual void multUTransInvToVector(bool withDiag, const Vector& in, Vector& out) const// (L [+D])^T)^-1^ * vector, нужно для бисопряженных градиентов
	{
		throw exception("У предобуславливателя хочу чтобы было multUTransInvToVector, а она есть только в одном классе");
	}

	virtual void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res) const = 0;  //Âåðõíèé
	virtual void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const = 0;
	virtual void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const = 0;     // (L [+D])^-1^ * vector
	virtual void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const = 0;     // (U [+D])^-1^ * vector
	virtual void getDiag(Vector& vec) const = 0; //Ïîëó÷èòü äèàãîíàëü
	// TODO Костыль от Сидорова А.
	virtual bool isSymmetric() const {
		return false;
	}
	virtual int Size() const = 0;
};

/*--------Автор базового класса, класса Вектор, и Разреженного Строчного формата - Ауст Кирилл----------------------*/

class MappedRowSkyLineFormat :public Matrix
{
	mutable vector<map<int, double>> supermatrix;
public:

	MappedRowSkyLineFormat(int size){
		supermatrix.resize(size);
	}
	MappedRowSkyLineFormat(bool t, istream& f) //Ввод матрицы, передаем имя файл, формат: n v(0,0)...v(n,n)
	{
		int n;
		double temp;
		f >> n;
		supermatrix.resize(n);
		if (t)
		{
			for (int i = 0; i < n; i++)
			{
				for (int j = 0; j < n; j++)
				{
					f >> temp;
					if (temp != 0)
					{
						supermatrix[i][j] = temp;
					}
				}
			}
		}
		else
		{
			int i, j;
			while (!f.eof())
			{
				f >> i >> j >> temp;
				SetElement(i, j, temp);
			}
		}
		if (!f)throw exception("End of matrix file read");
	}
	MappedRowSkyLineFormat(const Matrix& matr)
	{
		supermatrix.resize(matr.Size());
		matr.Run([&](int i, int j, double d){	SetElement(i, j, d); });
	}
	void Run(function<void(int, int, double)> fun) const override
	{
		for (unsigned int i = 0; i < supermatrix.size(); i++)
		{
			for (const auto& pair : supermatrix[i])
			{
				fun(i, pair.first, pair.second);
			}
		}
	}
	void MultTransMatrixVector(const Vector& vec, Vector& res) const override //Òðàíñïîíèðîâàííàÿ ìàòðèöà íà âåêòîð
	{
		res.resize(Size());
		for (auto& r : res)r = 0;
		Run([&](int i, int j, double d) {res[j] += vec[i] * d; });
	}
	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res) const override //Âåðõíèé òðåóãîëüíèê íà âåêòîð
	{// Òàê î÷åíü ìåäëåííî !!!!
		res.resize(Size());
		for(auto& r:res)r=0;

		Run([&](int i, int j, double d)
		{
			if (i == j && withDiag) res[i] += vec[i] * d;
			else
			{
				if (i<j) res[i] += vec[j] * d;
			}
		});
	}
	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override //Íèæíèé òðåóãîëüíèê íà âåêòîð
	{// Òàê î÷åíü ìåäëåííî !!!!
		res.resize(Size());
		for (auto& r : res)r = 0;

		Run([&](int i, int j, double d)
		{
			if (i == j && withDiag) res[i] += vec[i] * d;
			else
			{
				if (i>j) res[i] += vec[j] * d;
			}
		});
	}
	void SetElement(int i, int j, double d)
	{
		supermatrix[i][j] = d;
	}
	void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const override // (L [+D])^-1^ * vector
	{
		out.resize(Size());
		for (int i = 0; i < Size(); i++)
		{
			double temp = 0;
			for (const auto& pair : supermatrix[i])
			{
				if (pair.first < i)
				{
					temp += out[pair.first] * pair.second;
				}
			}
			if (withDiag) out[i] = (in[i] - temp) / supermatrix[i][i];
			else out[i] = (in[i] - temp);
		}
	}
	void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override // (U [+D])^-1^ * vector
	{
		out.resize(Size());
		for (int i = Size() - 1; i >= 0; i--)
		{
			double temp = 0;
			for (const auto& pair : supermatrix[i])
			{
				if (pair.first > i)
				{
					temp += out[pair.first] * pair.second;
				}
			}
			if (withDiag) out[i] = (in[i] - temp) / supermatrix[i][i];
			else out[i] = (in[i] - temp);
		}
	}


	virtual void multLTransInvToVector(bool withDiag, const Vector& in, Vector& out) const override// (L [+D])^T)^-1^ * vector
	{
		out.resize(Size());
		vector<double> temp;
		temp.resize(Size());

		for (int i = Size() - 1; i >= 0; i--)
		{
			out[i] = in[i] - temp[i];
			if (withDiag) out[i] /= supermatrix[i][i];

			for (const auto& pair : supermatrix[i])
			{
				if (pair.first < i)
				{
					temp[pair.first] += out[i] * pair.second;
				}
			}
		}
	}
	void multUTransInvToVector(bool withDiag, const Vector& in, Vector& out) const override// (U [+D])^T^-1^ * vector
	{
		out.resize(Size());
		vector<double> temp;
		temp.resize(Size());

		for (int i = 0; i < Size(); i++)
		{
			out[i] = in[i] - temp[i];
			if (withDiag) out[i] /= supermatrix[i][i];

			for (const auto& pair : supermatrix[i])
			{
				if (pair.first > i)
				{
					temp[pair.first] += out[i] * pair.second;
				}
			}
		}
	}

	void getDiag(Vector& vec) const override //Ïîëó÷èòü äèàãîíàëü
	{
		vec.resize(Size());
		for (unsigned int i = 0; i < supermatrix.size(); i++)
		{
			auto &temp = supermatrix[i].find(i);
			if (temp != supermatrix[i].end())
				vec[i] = temp->second;
			else vec[i] = 0;
		}
	}
	int Size() const override
	{
		return supermatrix.size();
	}

};

/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Профильный столбцово-строчный формат----------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Автор: Штро Роман-----------------------------------------------------------------------*/

class ProfileColumnRowFormat : public Matrix
{
private:
	vector<double> di, al, au;
	vector<int> ia;


public:
	ProfileColumnRowFormat(int size)
	{
		ia.resize(size + 1);
		di.resize(size);
	}
	ProfileColumnRowFormat(bool t, istream& f) //Ввод матрицы
	{
		int n, I, J;
		int k = 0;
		double temp;
		f >> n;
		vector<vector<double>> matr(n);

		if (t)
		{
			for (int i = 0; i < n; i++)
			{
				matr[i].resize(n);
				for (int j = 0; j < n; j++)
				{
					f >> temp;
					matr[i][j] = temp;
				}
			}
		}
		else
		{
			for (int i = 0; i < n; i++)
				matr[i].resize(n);
			while (!f.eof())
			{
				f >> I >> J >> temp;
				matr[I][J] = temp;
			}
		}

		ia.resize(n + 1);
		di.resize(n);
		ia[0] = 1;
		for (int i = 0; i < n; i++)
		{
			di[i] = matr[i][i];
			for (int j = n - 1; j > i && k == 0; j--)
			{
				if (matr[i][j] != 0 || matr[j][i] != 0)
				{
					k = 1;
					ia[i + 1] = ia[i] + j - i;
					for (int l = i + 1; l <= j; l++)
					{
						al.push_back(matr[l][i]);
						au.push_back(matr[i][l]);
					}
				}
			}
			if (k == 0)
				ia[i + 1] = ia[i];
			else
				k = 0;
		}
	}
	ProfileColumnRowFormat(const Matrix& matr)
	{
		int n;
		int k = 0;
		double temp;
		n = matr.Size();
		vector<vector<double>> matric(n);
		for (int i = 0; i < n; i++)
		{
			matric[i].resize(n);
		}
		matr.Run([&](int i, int j, double d){ matric[i][j] = d; });

		ia.resize(n + 1);
		di.resize(n);
		ia[0] = 1;
		for (int i = 0; i < n; i++)
		{
			di[i] = matric[i][i];
			for (int j = n - 1; j > i && k == 0; j--)
			{
				if (matric[i][j] != 0 || matric[j][i] != 0)
				{
					k = 1;
					ia[i + 1] = ia[i] + j - i;
					for (int l = i + 1; l <= j; l++)
					{
						al.push_back(matric[l][i]);
						au.push_back(matric[i][l]);
					}
				}
			}
			if (k == 0)
				ia[i + 1] = ia[i];
			else
				k = 0;
		}
	}
	void Run(function<void(int, int, double)> fun) const override
	{
		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (di[i] != 0)
				fun(i, i, di[i]);
			int k = i + 1;
			for (int j = 0; j < ia[i + 1] - ia[i]; j++)
			{
				if (al[l] != 0)
					fun(k + j, i, al[l]);
				if (au[l] != 0)
					fun(i, k + j, au[l]);
				l++;
			}
		}
	}
	/*void MultMatrixVector(const Vector& vec, Vector& res)
	{
	int l = 0;
	for (int i = 0; i < ia.size()-1; i++)
	{
	res[i] += di[i] * vec[i];
	int k = i + 1;
	for (int j = 0; j < ia[i + 1] - ia[i]; j++)
	{
	res[i] += au[l] * vec[k+j];
	res[k + j] += al[l] * vec[i];
	l++;
	}
	}
	}*/
	void MultTransMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			res[i] += di[i] * vec[i];
			int k = i + 1;
			for (int j = 0; j < ia[i + 1] - ia[i]; j++)
			{
				res[i] += al[l] * vec[k + j];
				res[k + j] += au[l] * vec[i];
				l++;
			}
		}
	}
	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res)const override   //Âåðõíèé
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (withDiag)
				res[i] += di[i] * vec[i];
			else
				res[i] += vec[i];
			int k = i + 1;
			for (int j = 0; j < ia[i + 1] - ia[i]; j++)
			{
				res[i] += au[l] * vec[k + j];
				l++;
			}
		}
	}
	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (withDiag)
				res[i] += di[i] * vec[i];
			else
				res[i] += vec[i];
			int k = i + 1;
			for (int j = 0; j < ia[i + 1] - ia[i]; j++)
			{
				res[k + j] += al[l] * vec[i];
				l++;
			}
		}
	}
	void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const override     // (L [+D])^-1^ * vector
	{
		out.resize(Size());
		out += in;
		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (withDiag)
				out[i] = out[i] / di[i];
			int k = i + 1;
			for (int j = 0; j < ia[i + 1] - ia[i]; j++)
			{
				out[k + j] -= al[l] * out[i];
				l++;
			}
		}
	}
	void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override     // (U [+D])^-1^ * vector
	{
		out.resize(Size());
		out += in;
		int l = au.size() - 1;
		for (int i = Size() - 1; i > 0; i--)
		{
			if (withDiag)
				out[i] = out[i] / di[i];
			for (int j = ia[i] - ia[i - 1] - 1; j >= 0; j--)
			{
				out[i - 1] -= au[l] * out[i + j];
				l--;
			}
		}
		if (withDiag)
			out[0] = out[0] / di[0];
	}
	void getDiag(Vector& vec) const override //Ïîëó÷èòü äèàãîíàëü
	{
		int n = di.size();
		vec.resize(n);
		for (int i = 0; i < n; i++)
			vec[i] = di[i];
	}
	int Size() const override
	{
		return di.size();
	}
};

/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Диагональный формат(произвольное количество диагоналей и отступов)----------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------Автор: Сарпова Юлия----------------------------------------------------------------------*/

class MappedDiagonalFormat :public Matrix
{
private:
	//размерность
	int n;
	//кол-во диагоналей
	int nd;
	//диагональ
	Vector di;
	//эл-ти верхнего треугольника
	vector <vector<double>> ggu;
	//эл-ти нижнего треугольника,
	vector <vector<double>> ggl;
	vector<int> ig;
	//вектор отступа - уменьшенные на 1 номера строк/столбцов ненулевых диагоналей
	map<int, int> mapig;//mapig[i] -- порядковый номер диагонали с номером i
	void Create(const Matrix& matr)
	{
		n = matr.Size();
		di.resize(n);
		//запись диагонали
		matr.getDiag(di);

		matr.Run([&](int i, int j, double d){
			if (i == j)return;
			int nd = abs(i - j);
			mapig[nd] = 0;
		});
		nd = 0;
		ig.resize(mapig.size());
		for (auto& el : mapig){ ig[nd] = el.first; el.second = nd++; }
		//запись вектор отступа диагоналей

		ggu.resize(nd);
		for (int i = 0; i < nd; i++)
			ggu[i].resize(n-ig[i]);

		ggl.resize(nd);
		for (int i = 0; i < nd; i++)
			ggl[i].resize(n - ig[i]);

		matr.Run([&](int i, int j, double d){
			if (i == j)return;
			if (j > i)
			{
				auto s = mapig[j - i];
				ggu[s][i] = d;
			}
			else
			{
				auto s = mapig[i - j];
				ggl[s][j] = d;
			}
		});
	}
public:

	MappedDiagonalFormat(const Matrix& matr)
	{
		Create(matr);
	}
	MappedDiagonalFormat(bool t, istream& f) //Ввод матрицы, передаем имя файл, формат: n v(0,0)...v(n,n)
	{
		auto input = make_shared<SpecialInputClass>(t, f);
		Create(*input);
	}
	// Из плотного формата в диагональный
	
	void Run(function<void(int, int, double)> fun)const override
	{
		int l = 0;
		for (int i = 0; i < n; i++)
		{
			if (di[i] != 0)
				fun(i, i, di[i]);
		}
		for (int i = 0; i < nd; i++)
			for (int j = 0; j < n - ig[i]; j++)
			{
			auto ir = j + ig[i];
			fun(ir, j, ggl[i][j]);
			fun(j, ir, ggu[i][j]);
			}
	}

	// Функция умножения матрицы на вектор, результат вектор
	void MultMatrixVector(const Vector& vec, Vector& res)const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int ir = 0;

		for (int i = 0; i < n; i++)
			res[i] = di[i] * vec[i];

		for (int i = 0; i < nd; i++)
		for (int j = 0; j < n - ig[i]; j++)
		{
			ir = j + ig[i];
			res[ir] += ggl[i][j] * vec[j];
			res[j] += ggu[i][j] * vec[ir];
		}

	}

	//Функция умножения транспонированной матрицы матрицы на вектор
	void MultTransMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int ir = 0;

		for (int i = 0; i < n; i++)
			res[i] = di[i] * vec[i];

		for (int i = 0; i < nd; i++)
		for (int j = 0; j < n - ig[i]; j++)
		{
			ir = j + ig[i];
			res[ir] += ggu[i][j] * vec[j];
			res[j] += ggl[i][j] * vec[ir];
		}

	}

	//функция умножения верхнего треугольника на вектор
	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res)  const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int ir = 0;

		for (int i = 0; i < n; i++)
			res[i] = 0;

		for (int i = 0; i < nd; i++)
		for (int j = 0; j < n - ig[i]; j++)
		{
			ir = j + ig[i];
			res[j] += ggu[i][j] * vec[ir];
		}
		if (withDiag)
		for (int i = 0; i < n; i++)
			res[i] += di[i] * vec[i];

	}

	//функция умножения нижнего треугольника на вектор
	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		int ir = 0;

		for (int i = 0; i < n; i++)
			res[i] = 0;

		for (int i = 0; i < nd; i++)
		for (int j = 0; j < n - ig[i]; j++)
		{
			ir = j + ig[i];
			res[ir] += ggl[i][j] * vec[j];
		}
		if (withDiag)
		for (int i = 0; i < n; i++)
			res[i] += di[i] * vec[i];
	}

	void multLInvToVector(bool withDiag, const Vector& in, Vector& out)  const override   // (L [+D])^-1^ * vector
	{

		out.resize(Size());

		for (int i = 1; i < n; i++)
			out[i] = 0;
		out[0] = in[0];
		for (int i = 1; i < n; i++)
		{
			for (int j = 0; j < i; j++)
			{
				auto find = mapig.find(i - j);
				if (find != mapig.end())
					out[i] += ggl[find->second][j] * out[j];
			}
			out[i] = (in[i] - out[i]);
			if (withDiag) out[i] = out[i] / di[i];
		}
	}
	void multUInvToVector(bool withDiag, const Vector& in, Vector& out)  const override   // (U [+D])^-1^ * vector
	{
		out.resize(Size());
		out += in;
		out[n - 1] = out[n - 1] / di[n - 1];
		for (int i = n - 2; i >= 0; i--)
		{
			for (int j = i + 1; j < n; j++)
			{
				auto find = mapig.find(j - i);
				if (find != mapig.end())
					out[i] -= ggu[find->second][i] * out[j];
			}
			out[i] = out[i] / di[i];
		}
	}

	//void SetElement(int i, int j, double d)
	//{
	//	auto top = [&](int i, int j, double d)
	//	{
	//		int s = 0;
	//		s = get_ig(j - i);
	//		ggu[s][i] = d;
	//	};

	//	auto lower = [&](int i, int j, double d)
	//	{
	//		int s = 0;
	//		s = get_ig(j - i);
	//		ggl[s][j] = d;
	//	};

	//	auto ddi = [&](int i, double d)
	//	{
	//		di[i] = d;
	//	};

	//	if (i < j) top(i, j, d);
	//	if (i > j) lower(i, j, d);
	//	if (i == j) ddi(i, d);
	//}

	void getDiag(Vector& vec)const override
	{
		vec = di;
	}
	int Size() const override
	{
		return n;
	}

};

/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Профильный Строчно-столбцовый формат----------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Автор: Нестеренко Тоня----------------------------------------------------------------------------------------*/

class ProfileRow_ColumnFormat : public Matrix
{
private:
	vector<double> di, al, au;//диагональ,нижний,верхний треугольник
	vector<int> ia; //массив профилЯ
public:
	ProfileRow_ColumnFormat(int size)
	{
		ia.resize(size + 1);
		di.resize(size);
	}

	ProfileRow_ColumnFormat(bool t, istream& f) //‚вод матрицы
	{
		int n, I, J;
		int k = 0;
		double temp;
		f >> n;
		vector<vector<double>> matr(n);
		//заполнЯем матрицу в плотный формат
		if (t)
		{
			for (int i = 0; i < n; i++)
			{
				matr[i].resize(n);
				for (int j = 0; j < n; j++)
				{
					f >> temp;
					matr[i][j] = temp;
				}
			}
		}
		else
		{
			for (int i = 0; i < n; i++)
				matr[i].resize(n);
			while (!f.eof())
			{
				f >> I >> J >> temp;
				matr[I][J] = temp;
			}
		}


		ia.resize(n + 1);
		di.resize(n);

		//заполнЯем диагональ и массив профилЯ ia
		for (int i = 0; i < n; i++)
			di[i] = matr[i][i];
		ia[0] = 1; ia[1] = 1; k = 0;
		for (int i = 1; i < n; i++)
		{
			for (int j = 0; j < i && k == 0; j++)
			{
				if (matr[i][j] != 0 || matr[j][i] != 0)
				{
					ia[i + 1] = ia[i] + (i - j);
					k = 1;
				}
				else { ia[i + 1] = ia[i]; }
			}
			k = 0;
		}
		int p = ia[n] - 1;
		al.resize(p);
		au.resize(p);

		//апись al и au
		for (int c = 1; c < n; c++)
		{
			int kol = ia[c + 1] - ia[c];
			if (kol != 0)
			{
				int j = c - kol;
				for (int i = ia[c]; i < ia[c + 1]; i++)
				{
					al[i - 1] = matr[c][j];
					au[i - 1] = matr[j][c];
					j++;
				}
			}
		}
	}



	//перевод из разреженного через плотный в профильный строчно-столбцовый
	ProfileRow_ColumnFormat(const Matrix& matr)
	{
		int n;
		int k = 0;
		double temp;
		n = matr.Size();
		vector<vector<double>> mat(n);
		for (int i = 0; i < n; i++)
		{
			mat[i].resize(n);
		}
		matr.Run([&](int i, int j, double d){ mat[i][j] = d; });

		ia.resize(n + 1);
		di.resize(n);

		//заполнЯем диагональ и массив профилЯ ia
		for (int i = 0; i < n; i++)
			di[i] = mat[i][i];
		ia[0] = 1; ia[1] = 1; k = 0;
		for (int i = 1; i < n; i++)
		{
			for (int j = 0; j < i && k == 0; j++)
			{
				if (mat[i][j] != 0 || mat[j][i] != 0)
				{
					ia[i + 1] = ia[i] + (i - j);
					k = 1;
				}
				else { ia[i + 1] = ia[i]; }
			}
			k = 0;
		}
		int p = ia[n] - 1;
		al.resize(p);
		au.resize(p);

		//апись al и au
		for (int c = 1; c < n; c++)
		{
			int kol = ia[c + 1] - ia[c];
			if (kol != 0)
			{
				int j = c - kol;
				for (int i = ia[c]; i < ia[c + 1]; i++)
				{
					al[i - 1] = mat[c][j];
					au[i - 1] = mat[j][c];
					j++;
				}
			}
		}
	}


	//пробег по ненулевым (длЯ разреж)
	void Run(function<void(int, int, double)> fun) const override
	{
		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (di[i] != 0)
				fun(i, i, di[i]);
			for (int j = i - (ia[i + 1] - ia[i]); j < i; j++)
			{
				if (al[l] != 0)
					fun(i, j, al[l]);
				if (au[l] != 0)
					fun(j, i, au[l]);
				l++;
			}

		}
	}
	//умножение транспонированной матрицы на вектор
	void MultTransMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int i = 0; i < ia.size() - 1; i++)
			res[i] += di[i] * vec[i];

		for (int i = 0; i < ia.size() - 1; i++)
		{
			int kol = ia[i + 1] - ia[i];
			if (kol != 0)
			{
				int nach = i - kol;
				int k = ia[i] - 1;
				for (int j = nach; j < i; j++)
				{
					res[i] = res[i] + au[k] * vec[j];
					res[j] = res[j] + al[k] * vec[i];
					k++;
				}
			}
		}
	}
	//умножение  матрицы на вектор
	void MultMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int i = 0; i < ia.size() - 1; i++)
			res[i] += di[i] * vec[i];

		for (int i = 0; i < ia.size() - 1; i++)
		{
			int kol = ia[i + 1] - ia[i];
			if (kol != 0)
			{
				int nach = i - kol;
				int k = ia[i] - 1;
				for (int j = nach; j < i; j++)
				{
					res[i] = res[i] + al[k] * vec[j];
					res[j] = res[j] + au[k] * vec[i];
					k++;
				}
			}
		}
	}
	//умножение нижнего на вектор 
	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (withDiag)res[i] = di[i] * vec[i];
			else  res[i] = vec[i];
		}

		for (int i = 0; i < ia.size() - 1; i++)
		{
			int kol = ia[i + 1] - ia[i];
			if (kol != 0)
			{
				int nach = i - kol;
				int k = ia[i] - 1;
				for (int j = nach; j < i; j++)
				{
					res[i] = res[i] + al[k] * vec[j];
					k++;
				}
			}
		}
	}
	//умножение верхнего треугольника на вектор
	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res)const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (withDiag)res[i] = di[i] * vec[i];
			else  res[i] = vec[i];
		}
		for (int i = 0; i < ia.size() - 1; i++)
		{
			int kol = ia[i + 1] - ia[i];
			if (kol != 0)
			{
				int nach = i - kol;
				int k = ia[i] - 1;
				for (int j = nach; j < i; j++)
				{
					res[j] = res[j] + au[k] * vec[i];
					k++;
				}
			}
		}
	}
	// (L [+D])^-1^ * vector
	void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const override
	{
		out.resize(Size());
		for (int i = 0; i < ia.size() - 1; i++)
		{
			double s = 0;
			int j = ia[i] - 1;
			int kol = ia[i + 1] - ia[i];

			for (int k = i - kol; k < i; k++)
			{
				s = s + out[k] * al[j];
				j++;
			}
			out[i] = (in[i] - s);
			if (withDiag) out[i] = out[i] / di[i];
		}
	}

	// (U [+D])^-1^ * vector
	void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override
	{
		out.resize(Size());
		vector<double> vs;
		vs.resize(Size());

		for (int i = 0; i < Size(); i++)
			vs[i] = in[i];

		for (int i = Size() - 1; i >= 0; i--)
		{
			double kk = 0;
			out[i] = vs[i];
			if (withDiag) out[i] = out[i] / di[i];
			int j = ia[i] - 1;
			int kol = ia[i + 1] - ia[i];

			for (int k = i - kol; k < i; k++)
			{
				kk = out[i] * au[j];
				vs[k] = vs[k] - kk;
				j++;
			}

		}

	}
	//получение диагонали
	void getDiag(Vector& vec) const override
	{
		int n = di.size();
		vec.resize(n);
		for (int i = 0; i < n; i++)
			vec[i] = di[i];
	}
	int Size() const override
	{
		return di.size();
	}
};


/*------------------------------Разреженный----------------------------*/
/*---------------------------столбцово-строчный------------------------*/
/*--------------------------------формат-------------------------------*/
class MappedRowColumnFormat : public Matrix
{
private:
	vector<double> di, al, au;//диагональ,нижний,верхний треугольник
	vector<int> ia, ja; //массив номеров стролбцов, индексов
	vector<vector<double>> matr;

public:
	MappedRowColumnFormat(int size)
	{
		ia.resize(size + 1);
		ja.resize(size + 1);
		di.resize(size);
	}

	MappedRowColumnFormat(bool t, istream& f) //ввод матрицы
	{
		int n, I, J;
		int k = 0;
		double temp;
		f >> n;
		matr.resize(n);


		if (t)
		{
			for (int i = 0; i < n; i++)
			{
				matr[i].resize(n);
				for (int j = 0; j < n; j++)
				{
					f >> temp;
					matr[i][j] = temp;
				}
			}
		}
		else
		{
			for (int i = 0; i < n; i++)
				matr[i].resize(n);
			while (!f.eof())
			{
				f >> I >> J >> temp;
				matr[I][J] = temp;
			}
		}


		ia.resize(n + 1);
		ja.resize(n*n);
		di.resize(n);


		//заполняем диагональ и массив профиля ia, ja
		for (int i = 0; i < n; i++)
			di[i] = matr[i][i];

		ia[0] = 1; ia[1] = 1; k = 0; int w = -1; int s = 1;
		for (int j = 0; j < n; j++)
		{
			for (int i = j + 1; i< n; i++)
			{
				//if(j!=i)
				//{
				if (matr[i][j] != 0 || matr[j][i] != 0)
				{
					w++;
					//ia[i + 1] = ia[i] + (i - j);
					//k = 1;
					s++;
					ia[j + 2] = s;
					ja[w] = i;

				}
				else { ia[i + 1] = ia[i]; }
				//}
			}
			//k = 0;
		}


		int p = ia[n] - 1;
		al.resize(p);
		au.resize(p);
		ja.resize(p);

		//запись al и au
		for (int c = 1; c < n; c++)
		{
			int kol = ia[c + 1] - ia[c];
			if (kol != 0)
			{

				//int j = abs(c - kol);
				for (int i = ia[c]; i < ia[c + 1]; i++)
				{
					int j = ja[i - 1];
					al[i - 1] = matr[j][c - 1];
					au[i - 1] = matr[c - 1][j];
					j++;
				}
			}
		}
	}



	//перевод через плотный в профильный строчно-столбцовый
	MappedRowColumnFormat(const Matrix& matr)
	{
		int n;
		int k = 0;
		double temp;
		n = matr.Size();
		vector<vector<double>> mat(n);
		for (int i = 0; i < n; i++)
		{
			mat[i].resize(n);
		}
		//matr.Run([&](int i, int j, double d){ mat[i][j] = d; });

		ia.resize(n + 1);
		ja.resize(n + 1);
		di.resize(n);

		//заполняем диагональ и массив профиля ia
		for (int i = 0; i < n; i++)
			di[i] = mat[i][i];
		ia[0] = 1; ia[1] = 1; k = 0;
		for (int i = 1; i < n; i++)
		{
			for (int j = 0; j < i && k == 0; j++)
			{
				if (mat[i][j] != 0 || mat[j][i] != 0)
				{
					ia[i + 1] = ia[i] + (i - j);
					k = 1;
				}
				else { ia[i + 1] = ia[i]; }
			}
			k = 0;
		}
		int p = ia[n] - 1;
		al.resize(p);
		au.resize(p);

		//запись al и au
		for (int c = 1; c < n; c++)
		{
			int kol = ia[c + 1] - ia[c];
			if (kol != 0)
			{
				int j = c - kol;
				for (int i = ia[c]; i < ia[c + 1]; i++)
				{
					al[i - 1] = mat[c][j];
					au[i - 1] = mat[j][c];
					j++;
				}
			}
		}
	}

	void Run(function<void(int, int, double)> fun) const override
	{
		int l = 0;
		for (int i = 0; i < ia.size() - 1; i++)
		{
			if (di[i] != 0)
				fun(i, i, di[i]);
			for (int j = i - (ia[i + 1] - ia[i]); j < i; j++)
			{
				if (al[l] != 0)
					fun(i, j, al[l]);
				if (au[l] != 0)
					fun(j, i, au[l]);
				l++;
			}

		}
	}
	//умножение транспонированной матрицы на вектор
	void MultTransMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int ii = 0; ii < ia.size() - 1; ii++)
			res[ii] = di[ii] * vec[ii];

		for (int j = 0; j < ia.size() - 1; j++)
		{
			int kol = ia[j + 1] - ia[j];
			if (kol != 0)
			{
				//int ii = i - kol;
				int kk = ia[j] - 1;
				int i = ja[j - 1] - 1;
				for (int l = 0; l < kol; l++)
				{
					//int i = ja[j-1];
					i++;
					res[i] = res[i] + au[kk] * vec[j - 1];
					res[j - 1] = res[j - 1] + al[kk] * vec[i];
					kk++;
				}
			}
		}
	}
	//умножение  матрицы на вектор
	void MultMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int ii = 0; ii < ia.size() - 1; ii++)
			res[ii] = di[ii] * vec[ii];

		for (int j = 0; j < ia.size() - 1; j++)
		{
			int kol = ia[j + 1] - ia[j];
			if (kol != 0)
			{
				//int ii = i - kol;
				int kk = ia[j] - 1;
				int i = ja[j - 1] - 1;
				for (int l = 0; l < kol; l++)
				{
					//int i = ja[j-1];
					i++;
					res[i] = res[i] + al[kk] * vec[j - 1];
					res[j - 1] = res[j - 1] + au[kk] * vec[i];
					kk++;
				}
			}
		}
	}
	//умножение нижнего треугольника на вектор 
	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int ii = 0; ii < ia.size() - 1; ii++)
			res[ii] = di[ii] * vec[ii];

		for (int j = 0; j < ia.size() - 1; j++)
		{
			int kol = ia[j + 1] - ia[j];
			if (kol != 0)
			{
				//int ii = i - kol;
				int kk = ia[j] - 1;
				int i = ja[j - 1] - 1;
				for (int l = 0; l < kol; l++)
				{
					//int i = ja[j-1];
					i++;
					res[i] = res[i] + al[kk] * vec[j - 1];
					//res[j-1] = res[j-1] + au[kk] * vec[i];
					kk++;
				}
			}
		}
	}
	//умножение верхнего треугольника на вектор
	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res)const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		for (int ii = 0; ii < ia.size() - 1; ii++)
			res[ii] = di[ii] * vec[ii];

		for (int j = 0; j < ia.size() - 1; j++)
		{
			int kol = ia[j + 1] - ia[j];
			if (kol != 0)
			{
				//int ii = i - kol;
				int kk = ia[j] - 1;
				int i = ja[j - 1] - 1;
				for (int l = 0; l < kol; l++)
				{
					//int i = ja[j-1];
					i++;
					//res[i] = res[i] + al[kk] * vec[j-1];
					res[j - 1] = res[j - 1] + au[kk] * vec[i];
					kk++;
				}
			}
		}
	}
	// (L [+D])^-1^ * vector
	void multLInvToVector(bool withDiag, const Vector& in, Vector& out)const override
	{
		out.resize(Size());
		for (int i = 0; i < ia.size() - 1; i++)
		{
			double s = 0;
			int j = ia[i] - 1;
			int kol = ia[i + 1] - ia[i];

			for (int k = abs(i - kol); k < i; k++)
			{
				s = s + out[k] * al[j];
				j++;
			}
			out[i] = (in[i] - s);
			if (withDiag) out[i] = out[i] / di[i];
		}
	}

	// (U [+D])^-1^ * vector
	/*void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override
	{
		out.resize(Size());
		vector<double> vs;
		vs.resize(Size());

		for (int i = 0; i < Size(); i++)
			vs[i] = in[i];

		for (int i = Size() - 1; i >= 0; i--)
		{
			double kk = 0;
			out[i] = vs[i];
			if (withDiag) out[i] = out[i] / di[i];
			int j = ia[i] - 1;
			int kol = ia[i + 1] - ia[i];

			for (int k = i - kol; k < i; k++)
			{
				kk = out[i] * au[j];
				vs[k] = vs[k] - kk;
				j++;
			}

		}

	}*/
	//получение диагонали
	void getDiag(Vector& vec) const override
	{
		int n = di.size();
		vec.resize(n);
		for (int i = 0; i < n; i++)
			vec[i] = di[i];
	}


	int Size() const override
	{
		return di.size();
	}
};
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Разреженный стобцовый формат-------------------------------------------------------------*/
/*-----------------------------------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------Автор: Пахоменко Дарья------------------------------------------------------------------*/

class MappedColumnSkyLineFormat :public Matrix
{
	mutable vector<map<int, double>> supermatrix;
public:

	MappedColumnSkyLineFormat(int size){
		supermatrix.resize(size);
	}

	MappedColumnSkyLineFormat(bool t, istream& f) //Ввод матрицы, передаем имя файл, формат: n v(0,0)...v(n,n)
	{
		int n;
		double temp;
		f >> n;
		supermatrix.resize(n);
		if (t) // ввод матрицы в плотном формате
		{
			for (int i = 0; i < n; i++)
			{
				for (int j = 0; j < n; j++)
				{
					f >> temp;
					if (temp != 0)
					{
						supermatrix[j][i] = temp;
					}
				}
			}
		}
		else // ввод элементов матрицы формата  i,j,elem
		{
			int i, j;
			while (!f.eof())
			{
				f >> i >> j >> temp;
				SetElement(i, j, temp);
			}
		}
		if (!f)throw exception("End of matrix file read");
	}

	MappedColumnSkyLineFormat(const Matrix& matr) // КОНВЕРТЕР в разреженный строчный формат
	{
		supermatrix.resize(matr.Size());
		matr.Run([&](int i, int j, double d){	SetElement(i, j, d); });
	}

	void Run(function<void(int, int, double)> fun) const override
	{
		for (unsigned int i = 0; i < supermatrix.size(); i++)
		{
			for (const auto& pair : supermatrix[i])
			{
				fun(pair.first, i, pair.second);
			}
		}
	}

	void MultTransMatrixVector(const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		Run([&](int j, int i, double d)
		{
			res[i] += vec[j] * d;
		});
	}

	void MultMatrixVector(const Vector& vec, Vector& res) const
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		Run([&](int j, int i, double d)
		{
			res[j] += vec[i] * d;
		});
	}

	void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		Run([&](int i, int j, double d)
		{
			if (i == j && withDiag) res[i] += vec[i] * d;
			else
			{
				if (i>j)
					res[i] += vec[j] * d;
			}
		});
	}

	void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res) const override
	{
		res.resize(Size());
		for (auto& r : res)r = 0;

		Run([&](int i, int j, double d)
		{
			if (i == j && withDiag) res[i] += vec[i] * d;
			else
			{
				if (i<j)
					res[i] += vec[j] * d;
			}
		});
	}

	void SetElement(int i, int j, double d)
	{
		supermatrix[j][i] = d;
	}

	void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override // (U [+D])^-1^ * vector
	{
		out.resize(Size());
		vector<double> temp;
		temp.resize(Size());

		for (int i = Size() - 1; i >= 0; i--)
		{
			out[i] = in[i] - temp[i];
			if (withDiag) out[i] /= supermatrix[i][i];

			for (const auto& pair : supermatrix[i])
			{
				if (pair.first < i)
				{
					temp[pair.first] += out[i] * pair.second;
				}
			}
		}
	}

	void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const override // (L [+D])^-1^ * vector
	{
		out.resize(Size());
		vector<double> temp;
		temp.resize(Size());

		for (int i = 0; i < Size(); i++)
		{
			out[i] = in[i] - temp[i];
			if (withDiag) out[i] /= supermatrix[i][i];

			for (const auto& pair : supermatrix[i])
			{
				if (pair.first > i)
				{
					temp[pair.first] += out[i] * pair.second;
				}
			}
		}
	}

	void getDiag(Vector& vec) const override
	{
		vec.resize(Size());
		for (unsigned int i = 0; i < supermatrix.size(); i++)
		{
			auto &temp = supermatrix[i].find(i);
			if (temp != supermatrix[i].end())
				vec[i] = temp->second;
			else vec[i] = 0;
		}
	}

	int Size() const override
	{
		return supermatrix.size();
	}

};